%% mf_analyze, print and plot results fo simulation

%% Analyzing output of xsecton model
% Plot /shows the results of the smulation for a given stress period and
% time step to be set by the user
% TO 090223 111212
% Bas des Tombe

clear variables
close all

%% load model name and basename contained in name.mat
load name
load(basename);  % this yields the stored matrices that make up the grid and the model
load underneath
load setup
load results
load CurrentSeriesNumber.mat
n=CurNum;  % Current setup number

%% Start
[NRow,NCol,NLay]=size(IBOUND);

%% compute coordinates of axes and cell centers
[xGr,yGr,zGr,xm,ym,zm,Dx,Dy,Dz,Nx,Ny,Nz]=modelsize3(xGr,yGr,zGr);

%% load the unformatted files with the heads, the concentration and the budget terms
H=readDat([basename,'.hds']);
H=maskHC(H,IBOUND & H(1).values<1000);   % read heads and mask inactive cells with NaN

C=readMT3D('MT3D001.UCN');
C=maskHC(C,IBOUND);   % read concentrations and mask inactive with NaN
C=maskHC(C,ICBUND);   % read concentrations and mask inactive with NaN

B=readBud([basename,'.bgt'],'FLOWRIGHTFACE');  B=maskHC(B,IBOUND); % read budgetfile and mask

for i=1:length(B)
  for j=1:length(B(i).term)
      B(i).term{j}(IBOUND==0 | H(1).values>1000)=NaN;
  end
end

B=mf_Psi(B);
dPsi=1;
crange=ContourRange(B,dPsi,'','Psi');
minPsi=min(crange);
maxPsi=max(crange);

% Stream function mask, to not plot stream lines above sand body
PSIMASK=NaN(Nz+1,Nx-1);
PSIMASK(inpolygon(ones(size(zGr(:)))*xGr(2:end-1),zGr(:)*ones(size(xGr(2:end-1))),xSand,zSand))=1;
% PSIMASK(inpolygon(ones(size(zGr(:)))*xGr(2:end-1),zGr(:)*ones(size(xGr(2:end-1))),xCanL,zCanL))=NaN;

concs= -0.05:0.1:1.05;      % Conc is always between 0 and 1
chead=0:0.05:65;      % Heads is always between base and top of model
cpsi =minPsi:(maxPsi-minPsi)/50:maxPsi; % Psi-range is computed above

%% Freshwater stored
Quantity=[0;0];
ConcentrationWell=[0;0];
ConcentrationMonitorPoint1=[0;0];
MaxQuantity=dy*dx*ixWall*dz*izWall*peff;
MaxAllowedC=150;    % [mg/l]
ExtractionStart=C(find(QMilk<0,1,'first')-1).time;    % Time start extracting[min]
CheckMonitorPoint1=1;
CheckMonitorPoint2=1;
dummy=0;
MaxUsableCol=0;
SizeTooSmall=0;
ConvertGramsSalt=29000; % Convert relative concentration to milligrams of salt per liter


for i=1:length(C);
    A=XS(C(i).values).*ConvertGramsSalt; % convert concentration values
    Quantity=[Quantity, [XS(C(i).time);length(find(A(1:izWall,1:ixWall)<(MaxAllowedC)))*dy*dx*dz*peff]];    % [cm3]

%%  Concentration Monitor Point 1
    if CheckMonitorPoint1        
        ConcentrationMonitorPoint1=[ConcentrationMonitorPoint1,[XS(C(i).time);A(izWall, ixWall)]];   % Concentration at Well
        if ConcentrationMonitorPoint1(end,end)<=0.95*ConvertGramsSalt;
            SizeTooSmall=1;
            SizeTooSmallTime=ConcentrationMonitorPoint1(1,end);
            CheckMonitorPoint1=0;
        end
    end

%%  Concentration Well
    ConcentrationWell=[ConcentrationWell, [XS(C(i).time);A(izMilk, ixMilk)]]; % Concentration at Well
    if (ConcentrationWell(end,end)>=MaxAllowedC && XS(C(i).time)>=ExtractionStart && dummy==0); % Log when extracting is not allowed anymore, exceeds MaxAllowedC, just ones
        APrev=XS(C(i-1).values); % rename previous concentration values
        MaxUsableQuantity=[XS(C(i-1).time);length(find(APrev(1:izWall,1:ixWall)<MaxAllowedC))*dy*dx*dz];
        CurConcentration=ConcentrationWell(:,end);
        MaxUsableCol=i+1;
        CheckMonitorPoint2=0;
        dummy=1;
    end
end

%% Recovery efficiency
Injection=[Quantity(1,:);0,QMilk'.*diff(ConcentrationWell(1,:))*MULTIMILK(n)];  % Injection[cm3/delta t],QMilk[l/h],Diffcon[min]

[m]=find(cumsum(Injection(2,:))==max(cumsum(Injection(2,:))),1,'first');  %Col where injection stoppes
[o]=find(cumsum(Injection(2,:))==max(cumsum(Injection(2,:))),1,'last');  %Col where dejection startes
QInjected=[cumsum(Injection(2,1:m)), ones(1,length(C)+1-m)*max(cumsum(Injection(2,1:m)))];
QRetreived=[zeros(1,o),cumsum(Injection(2,o:length(C)))];

% if (ExtractionStart>=ConcentrationWell(1,o) || m==0 || o==0)
%     error('Please set ExtractionStart time=%d[min], in mf_analyzepretty on line 84.\n This error should be removed',ConcentrationWell(1,m));
% end

RE=(QRetreived(MaxUsableCol)/QInjected(MaxUsableCol))*-100;

fprintf('Recovery efficiency is expected to be %2.1f percent \n',RE);


%% Save as file
ConcentrationSaved=XS(C(MaxUsableCol).values);
save ConcSaved ConcentrationSaved;    
fprintf('Concentration at stop point saved in "ConcSaved.mat".\n')

if SizeTooSmall,
    SizeTooSmallSaved(n,1)=SizeTooSmallTime;
else
    SizeTooSmallSaved(n,1)=0;
end

ConcentrationWellSaved(n,1)=ConcentrationWell(1,m);  % After t=x[min], start retreiving
QInjectedSaved(n,1)=QInjected(m);  % Total injected
QRetreivedSaved(n,1)=(-1*QRetreived(1,MaxUsableCol(1)));  % Total retreived
RESaved(n,1)=RE;  % Max recovery, untill conc exceedes maxconc
MaxQuantitySaved(n,1)=MaxQuantity;  % Max storeable freshwater

save results ConcentrationWellSaved QInjectedSaved QRetreivedSaved RESaved SizeTooSmallSaved MaxQuantity MaxQuantitySaved -append;